source('misc_functions/functions.R')
library(tidyverse)
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.8
✔ tidyr   0.8.2     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(plotly)

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout
library(modelr)
library(mgcv)
Loading required package: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:dplyr’:

    collapse

This is mgcv 1.8-25. For overview type 'help("mgcv-package")'.

The data set has been constructed using average Science scores by country from the Programme for International Student Assessment (PISA) 2006, along with GNI per capita (Purchasing Power Parity, 2005 dollars), Educational Index, Health Index, and Human Development Index from UN data. The key variables are as follows:

The first thing to do is get the data in and do some initial inspections.

We can look at the individual relationships of covariates with overall science score.

dmelt = pisa %>% 
  select(-Evidence, -Explain, -Issues) %>% 
  gather(key=Variable, 
         value=Value, 
         -Overall, -Country)

ggplot(aes(x=Value,y=Overall), data=dmelt) +
  geom_point(color='#ff5500',alpha=.75) +
  geom_smooth(se=F, lwd=.5, color='#00aaff') +
  geom_text(aes(label=Country), alpha=0, size=1,angle=30, hjust=-.2,   # making transparent so only plotly will show the country
            vjust=-.2) +
  facet_wrap(~Variable, scales='free_x') +
  labs(x='') +
  theme_trueMinimal()

Single Predictor

Linear Fit

We will start with the simple situation of a single predictor. Let’s begin by using a typical linear regression to predict science scores by the Income index. We could use the standard R lm function, but I’ll leave that as an exercise for comparison. We can still do straightforward linear models with the gam function, and again it is important to note that the standard linear model can be seen as a special case of a GAM.


Family: gaussian 
Link function: identity 

Formula:
Overall ~ Income

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   204.32      35.37   5.777 4.32e-07 ***
Income        355.85      46.79   7.606 5.36e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


R-sq.(adj) =  0.518   Deviance explained = 52.7%
GCV = 1504.5  Scale est. = 1448.8    n = 54

GAM

Fitting the model

Let’s now fit an actual generalized additive model using the same cubic spline as our basis. We again use the gam function as before for basic model fitting, but now we are using a function s within the formula to denote the smooth terms. Within that function we also specify the type of smooth, though a default is available. I chose bs = cr, denoting cubic regression splines.

mod_gam1 <- gam(Overall ~ s(Income, bs="cr"), data=pisa)
summary(mod_gam1)

Family: gaussian 
Link function: identity 

Formula:
Overall ~ s(Income, bs = "cr")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  470.444      4.082   115.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df     F  p-value    
s(Income) 6.895  7.741 16.67 1.59e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =    0.7   Deviance explained = 73.9%
GCV = 1053.7  Scale est. = 899.67    n = 54

In the summary, we first see the distribution assumed, as well as the link function used, in this case normal and identity, respectively, which to iterate, had we had no smoothing, would result in a standard linear model.

We have two general pieces of output:

  • Parametric: any non-smooth term, interpreted as any coefficient in the same GLM model.
  • Smooth terms: stuff we allow to ‘wiggle’ (or otherwise have a special relationship), interpreted visually via component plot.

In this case, the only parametric component is the intercept, but it’s good to remember that you are not bound to smooth every effect of interest, and indeed, as we will discuss in more detail later, part of the process may involve refitting the model with terms that were found to be linear for the most part anyway. The smooth component of our model regarding a country’s income and its relationship with overall science score suggests that it is statistically significant, but there are a couple of things in the model summary that would be unfamiliar. You may consult the primary document for details, but the main thing is that you can interpret it like any result you’d see in a regression table. The statistical significance is an approximation though, so again, let visualization be your guide.

Since this is in the end a linear model with a normal distribution for the target, we can use the R^2 value as we would with lm. Note that it already is adjusted for you. The deviance explained is the unadjusted R^2 in this case, but generalizes to other models that don’t assume the normal distribution. The scale estimate is our residual sums of squares from the standard regression setting, and the GCV is a better estimate of what that would be on held-out (i.e. new) data. The point is that you are seeing the usual output you have with other models, just with different labels.

Graphical Display

The primary way to interpret the effect of income is visually though. The mgcv package will provide a basic plot as follows:

plot(mod_gam1)

What you’re seeing is a component plot with the y/fitted values centered. Almost no one finds this intuitive in practice, nor is it usually what they’d report. Instead we can get predictions at key values of other covariates (e.g. held at their mean or reference category). In this case with only one predictor, we just get the predicted values on the original scale of the outcome. We can use the ggeffects package for this (or see my visibly package).

library(ggeffects)

plot_dat <- ggpredict(mod_gam1, terms = "Income")

ggplot(plot_dat, aes(x = x, y = predicted)) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(color = 'dodgerblue') + 
  labs(x = 'Income')

Model Comparison

We can compare models via AIC or the GCV coefficient, both of which focus on the predictive capabilities of the model on new data.

AIC(mod_lm, mod_gam1)

Likelihood ratio test (approximate).

anova(mod_lm, mod_gam1, test="Chisq")
Analysis of Deviance Table

Model 1: Overall ~ Income
Model 2: Overall ~ s(Income, bs = "cr")
  Resid. Df Resid. Dev     Df Deviance  Pr(>Chi)    
1    52.000      75336                              
2    45.259      41479 6.7411    33857 2.778e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple Predictors

Linear FIt

mod_lm2 <- gam(Overall ~ Income + Edu + Health, data=pisa)
summary(mod_lm2)

Family: gaussian 
Link function: identity 

Formula:
Overall ~ Income + Edu + Health

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   121.18      78.97   1.535   0.1314    
Income        182.32      85.27   2.138   0.0376 *  
Edu           234.11      54.78   4.274 9.06e-05 ***
Health         27.01     134.90   0.200   0.8421    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


R-sq.(adj) =  0.616   Deviance explained = 63.9%
GCV = 1212.3  Scale est. = 1119      n = 52

GAM

mod_gam2 <- gam(Overall ~ s(Income) + s(Edu) + s(Health), data=pisa)
summary(mod_gam2)

Family: gaussian 
Link function: identity 

Formula:
Overall ~ s(Income) + s(Edu) + s(Health)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  471.154      2.772     170   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df     F  p-value    
s(Income) 7.593  8.415 8.826 1.33e-07 ***
s(Edu)    6.204  7.178 3.308  0.00733 ** 
s(Health) 1.000  1.000 2.736  0.10661    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.863   Deviance explained = 90.3%
GCV = 573.83  Scale est. = 399.5     n = 52

We can see the fit appears to be better. In addition, as the effective degrees of freedom is 1 for Health, it is essentially a linear fit, so if desired, we can leave it among the parametric terms.

# plot(mod_gam2)

library(patchwork) # devtools::install_github("thomasp85/patchwork")

g1 = 
  ggpredict(mod_gam2, terms = "Income") %>% 
  ggplot(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(color = 'dodgerblue') + 
  labs(x = 'Income')
g2 = 
  ggpredict(mod_gam2, terms = "Edu") %>% 
  ggplot(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(color = 'dodgerblue') + 
  labs(x = 'Edu')
g3 = 
  ggpredict(mod_gam2, terms = "Health") %>% 
  ggplot(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  geom_line(color = 'dodgerblue') + 
  labs(x = 'Health')

(g2 + g3 + g1 + plot_layout(nrow = 2)) * theme_trueMinimal()

2D Smooths

mod_gam3 <- gam(Overall ~ Health + te(Income, Edu), data=pisa)
summary(mod_gam3)

Family: gaussian 
Link function: identity 

Formula:
Overall ~ Health + te(Income, Edu)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    574.5      101.4   5.666 1.42e-06 ***
Health        -115.8      113.5  -1.020    0.314    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                 edf Ref.df     F  p-value    
te(Income,Edu) 10.31  12.41 9.084 2.43e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.804   Deviance explained = 84.8%
GCV = 747.52  Scale est. = 570.59    n = 52
# requires the visibly package, otherwise see the main document; 
# devtools::install_github("m-clark/visibly")
visibly::plot_gam_3d(model = mod_gam3, main_var = Income, second_var = Edu, palette='bilbao', direction = 1)

Model Comparison

AIC(mod_lm2, mod_gam2)
LS0tCnRpdGxlOiAiQXBwbGljYXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawplZGl0b3Jfb3B0aW9uczogCiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQotLS0KCmBgYHtyIG1pc2NfZnVuY3Rpb25zfQpzb3VyY2UoJ21pc2NfZnVuY3Rpb25zL2Z1bmN0aW9ucy5SJykKYGBgCgpgYGB7ciBsb2FkX3BhY2thZ2VzfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShwbG90bHkpCmxpYnJhcnkobW9kZWxyKQpsaWJyYXJ5KG1nY3YpCmBgYAoKVGhlIGRhdGEgc2V0IGhhcyBiZWVuIGNvbnN0cnVjdGVkIHVzaW5nIGF2ZXJhZ2UgU2NpZW5jZSBzY29yZXMgYnkgY291bnRyeSBmcm9tIHRoZSBQcm9ncmFtbWUgZm9yIEludGVybmF0aW9uYWwgU3R1ZGVudCBBc3Nlc3NtZW50IChQSVNBKSAyMDA2LCBhbG9uZyB3aXRoIEdOSSBwZXIgY2FwaXRhIChQdXJjaGFzaW5nIFBvd2VyIFBhcml0eSwgMjAwNSBkb2xsYXJzKSwgRWR1Y2F0aW9uYWwgSW5kZXgsIEhlYWx0aCBJbmRleCwgYW5kIEh1bWFuIERldmVsb3BtZW50IEluZGV4IGZyb20gVU4gZGF0YS4gVGhlIGtleSB2YXJpYWJsZXMgYXJlIGFzIGZvbGxvd3M6CgoKLSBPdmVyYWxsIFNjaWVuY2UgU2NvcmUgKGF2ZXJhZ2Ugc2NvcmUgZm9yIDE1IHllYXIgb2xkcykKLSBJbnRlcmVzdCBpbiBzY2llbmNlCi0gU3VwcG9ydCBmb3Igc2NpZW50aWZpYyBpbnF1aXJ5Ci0gSW5jb21lIEluZGV4Ci0gSGVhbHRoIEluZGV4Ci0gRWR1Y2F0aW9uIEluZGV4Ci0gSHVtYW4gRGV2ZWxvcG1lbnQgSW5kZXggKGNvbXBvc2VkIG9mIHRoZSBJbmNvbWUgaW5kZXgsIEhlYWx0aCBJbmRleCwgYW5kIEVkdWNhdGlvbiBJbmRleCkKClRoZSBmaXJzdCB0aGluZyB0byBkbyBpcyBnZXQgdGhlIGRhdGEgaW4gYW5kIGRvIHNvbWUgaW5pdGlhbCBpbnNwZWN0aW9ucy4KCmBgYHtyIGluaXRpYWxfaW5zcGVjdGlvbl9vZl9waXNhLCBlY2hvPTF9CnBpc2EgPSByZWFkLmNzdignZGF0YS9waXNhc2NpMjAwNi5jc3YnKQoKcGlzYSAKCnBpc2EgJT4lIHNlbGVjdCgtQ291bnRyeSkKYGBgCgpXZSBjYW4gbG9vayBhdCB0aGUgaW5kaXZpZHVhbCByZWxhdGlvbnNoaXBzIG9mIGNvdmFyaWF0ZXMgd2l0aCBvdmVyYWxsIHNjaWVuY2Ugc2NvcmUuCgpgYGB7ciBiaXZhcmlhdGVfcmVsYXRpb25zaGlwcywgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KZG1lbHQgPSBwaXNhICU+JSAKICBzZWxlY3QoLUV2aWRlbmNlLCAtRXhwbGFpbiwgLUlzc3VlcykgJT4lIAogIGdhdGhlcihrZXk9VmFyaWFibGUsIAogICAgICAgICB2YWx1ZT1WYWx1ZSwgCiAgICAgICAgIC1PdmVyYWxsLCAtQ291bnRyeSkKCmdncGxvdChhZXMoeD1WYWx1ZSx5PU92ZXJhbGwpLCBkYXRhPWRtZWx0KSArCiAgZ2VvbV9wb2ludChjb2xvcj0nI2ZmNTUwMCcsYWxwaGE9Ljc1KSArCiAgZ2VvbV9zbW9vdGgoc2U9RiwgbHdkPS41LCBjb2xvcj0nIzAwYWFmZicpICsKICBnZW9tX3RleHQoYWVzKGxhYmVsPUNvdW50cnkpLCBhbHBoYT0wLCBzaXplPTEsYW5nbGU9MzAsIGhqdXN0PS0uMiwgICAjIG1ha2luZyB0cmFuc3BhcmVudCBzbyBvbmx5IHBsb3RseSB3aWxsIHNob3cgdGhlIGNvdW50cnkKICAgICAgICAgICAgdmp1c3Q9LS4yKSArCiAgZmFjZXRfd3JhcCh+VmFyaWFibGUsIHNjYWxlcz0nZnJlZV94JykgKwogIGxhYnMoeD0nJykgKwogIHRoZW1lX3RydWVNaW5pbWFsKCkKYGBgCgoKIyMgU2luZ2xlIFByZWRpY3RvcgoKIyMjIExpbmVhciBGaXQKCldlIHdpbGwgc3RhcnQgd2l0aCB0aGUgc2ltcGxlIHNpdHVhdGlvbiBvZiBhIHNpbmdsZSBwcmVkaWN0b3IuIExldCdzIGJlZ2luIGJ5IHVzaW5nIGEgdHlwaWNhbCBsaW5lYXIgcmVncmVzc2lvbiB0byBwcmVkaWN0IHNjaWVuY2Ugc2NvcmVzIGJ5IHRoZSBJbmNvbWUgaW5kZXguICBXZSBjb3VsZCB1c2UgdGhlIHN0YW5kYXJkIFIgYGxtYCBmdW5jdGlvbiwgYnV0IEknbGwgbGVhdmUgdGhhdCBhcyBhbiBleGVyY2lzZSBmb3IgY29tcGFyaXNvbi4gIFdlIGNhbiBzdGlsbCBkbyBzdHJhaWdodGZvcndhcmQgbGluZWFyIG1vZGVscyB3aXRoIHRoZSBgZ2FtYCBmdW5jdGlvbiwgYW5kIGFnYWluIGl0IGlzIGltcG9ydGFudCB0byBub3RlIHRoYXQgdGhlIHN0YW5kYXJkIGxpbmVhciBtb2RlbCBjYW4gYmUgc2VlbiBhcyBhIHNwZWNpYWwgY2FzZSBvZiBhIEdBTS4KCmBgYHtyIG1vZF9sbSwgZWNobz0tNH0KbGlicmFyeShtZ2N2KQptb2RfbG0gPC0gZ2FtKE92ZXJhbGwgfiBJbmNvbWUsIGRhdGE9cGlzYSkKc3VtbWFyeShtb2RfbG0pCmBgYAoKIyMjIEdBTQoKIyMjIEZpdHRpbmcgdGhlIG1vZGVsCgpMZXQncyBub3cgZml0IGFuIGFjdHVhbCBnZW5lcmFsaXplZCBhZGRpdGl2ZSBtb2RlbCB1c2luZyB0aGUgc2FtZSBjdWJpYyBzcGxpbmUgYXMgb3VyIGJhc2lzLiBXZSBhZ2FpbiB1c2UgdGhlIGBnYW1gIGZ1bmN0aW9uIGFzIGJlZm9yZSBmb3IgYmFzaWMgbW9kZWwgZml0dGluZywgYnV0IG5vdyB3ZSBhcmUgdXNpbmcgYSBmdW5jdGlvbiBgc2Agd2l0aGluIHRoZSBmb3JtdWxhIHRvIGRlbm90ZSB0aGUgc21vb3RoIHRlcm1zLiAgV2l0aGluIHRoYXQgZnVuY3Rpb24gd2UgYWxzbyBzcGVjaWZ5IHRoZSB0eXBlIG9mIHNtb290aCwgdGhvdWdoIGEgZGVmYXVsdCBpcyBhdmFpbGFibGUuICBJIGNob3NlIGBicyA9IGNyYCwgZGVub3RpbmcgY3ViaWMgcmVncmVzc2lvbiBzcGxpbmVzLgoKYGBge3IgbW9kX2dhbTF9Cm1vZF9nYW0xIDwtIGdhbShPdmVyYWxsIH4gcyhJbmNvbWUsIGJzPSJjciIpLCBkYXRhPXBpc2EpCnN1bW1hcnkobW9kX2dhbTEpCmBgYAoKSW4gdGhlIHN1bW1hcnksIHdlIGZpcnN0IHNlZSB0aGUgZGlzdHJpYnV0aW9uIGFzc3VtZWQsIGFzIHdlbGwgYXMgdGhlIGxpbmsgZnVuY3Rpb24gdXNlZCwgaW4gdGhpcyBjYXNlIG5vcm1hbCBhbmQgaWRlbnRpdHksIHJlc3BlY3RpdmVseSwgd2hpY2ggdG8gaXRlcmF0ZSwgaGFkIHdlIGhhZCBubyBzbW9vdGhpbmcsIHdvdWxkIHJlc3VsdCBpbiBhIHN0YW5kYXJkIGxpbmVhciBtb2RlbC4KCldlIGhhdmUgdHdvIGdlbmVyYWwgcGllY2VzIG9mIG91dHB1dDoKCi0gKipQYXJhbWV0cmljKio6IGFueSBub24tc21vb3RoIHRlcm0sIGludGVycHJldGVkIGFzIGFueSBjb2VmZmljaWVudCBpbiB0aGUgc2FtZSBHTE0gbW9kZWwuCi0gKipTbW9vdGggdGVybXMqKjogc3R1ZmYgd2UgYWxsb3cgdG8gJ3dpZ2dsZScgKG9yIG90aGVyd2lzZSBoYXZlIGEgc3BlY2lhbCByZWxhdGlvbnNoaXApLCBpbnRlcnByZXRlZCB2aXN1YWxseSB2aWEgY29tcG9uZW50IHBsb3QuCgpJbiB0aGlzIGNhc2UsIHRoZSBvbmx5IHBhcmFtZXRyaWMgY29tcG9uZW50IGlzIHRoZSBpbnRlcmNlcHQsIGJ1dCBpdCdzIGdvb2QgdG8gcmVtZW1iZXIgdGhhdCB5b3UgYXJlIG5vdCBib3VuZCB0byBzbW9vdGggZXZlcnkgZWZmZWN0IG9mIGludGVyZXN0LCBhbmQgaW5kZWVkLCBhcyB3ZSB3aWxsIGRpc2N1c3MgaW4gbW9yZSBkZXRhaWwgbGF0ZXIsIHBhcnQgb2YgdGhlIHByb2Nlc3MgbWF5IGludm9sdmUgcmVmaXR0aW5nIHRoZSBtb2RlbCB3aXRoIHRlcm1zIHRoYXQgd2VyZSBmb3VuZCB0byBiZSBsaW5lYXIgZm9yIHRoZSBtb3N0IHBhcnQgYW55d2F5LiAgVGhlIHNtb290aCBjb21wb25lbnQgb2Ygb3VyIG1vZGVsIHJlZ2FyZGluZyBhIGNvdW50cnkncyBpbmNvbWUgYW5kIGl0cyByZWxhdGlvbnNoaXAgd2l0aCBvdmVyYWxsIHNjaWVuY2Ugc2NvcmUgc3VnZ2VzdHMgdGhhdCBpdCBpcyBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50LCBidXQgdGhlcmUgYXJlIGEgY291cGxlIG9mIHRoaW5ncyBpbiB0aGUgbW9kZWwgc3VtbWFyeSB0aGF0IHdvdWxkIGJlIHVuZmFtaWxpYXIuICBZb3UgbWF5IGNvbnN1bHQgdGhlIHByaW1hcnkgZG9jdW1lbnQgZm9yIGRldGFpbHMsIGJ1dCB0aGUgbWFpbiB0aGluZyBpcyB0aGF0IHlvdSBjYW4gaW50ZXJwcmV0IGl0IGxpa2UgYW55IHJlc3VsdCB5b3UnZCBzZWUgaW4gYSByZWdyZXNzaW9uIHRhYmxlLiAgVGhlIHN0YXRpc3RpY2FsIHNpZ25pZmljYW5jZSBpcyBhbiBhcHByb3hpbWF0aW9uIHRob3VnaCwgc28gYWdhaW4sIGxldCB2aXN1YWxpemF0aW9uIGJlIHlvdXIgZ3VpZGUuCgpTaW5jZSB0aGlzIGlzIGluIHRoZSBlbmQgYSBsaW5lYXIgbW9kZWwgd2l0aCBhIG5vcm1hbCBkaXN0cmlidXRpb24gZm9yIHRoZSB0YXJnZXQsIHdlIGNhbiB1c2UgdGhlIFJeMiB2YWx1ZSBhcyB3ZSB3b3VsZCB3aXRoIGBsbWAuICBOb3RlIHRoYXQgaXQgYWxyZWFkeSBpcyBhZGp1c3RlZCBmb3IgeW91LiBUaGUgKmRldmlhbmNlIGV4cGxhaW5lZCogaXMgdGhlIHVuYWRqdXN0ZWQgUl4yIGluIHRoaXMgY2FzZSwgYnV0IGdlbmVyYWxpemVzIHRvIG90aGVyIG1vZGVscyB0aGF0IGRvbid0IGFzc3VtZSB0aGUgbm9ybWFsIGRpc3RyaWJ1dGlvbi4gVGhlICpzY2FsZSBlc3RpbWF0ZSogaXMgb3VyIHJlc2lkdWFsIHN1bXMgb2Ygc3F1YXJlcyBmcm9tIHRoZSBzdGFuZGFyZCByZWdyZXNzaW9uIHNldHRpbmcsIGFuZCB0aGUgKkdDViogaXMgYSBiZXR0ZXIgZXN0aW1hdGUgb2Ygd2hhdCB0aGF0IHdvdWxkIGJlIG9uIGhlbGQtb3V0IChpLmUuIG5ldykgZGF0YS4gIFRoZSBwb2ludCBpcyB0aGF0IHlvdSBhcmUgc2VlaW5nIHRoZSB1c3VhbCBvdXRwdXQgeW91IGhhdmUgd2l0aCBvdGhlciBtb2RlbHMsIGp1c3Qgd2l0aCBkaWZmZXJlbnQgbGFiZWxzLgoKIyMjIEdyYXBoaWNhbCBEaXNwbGF5CgpUaGUgcHJpbWFyeSB3YXkgdG8gaW50ZXJwcmV0IHRoZSBlZmZlY3Qgb2YgaW5jb21lIGlzIHZpc3VhbGx5IHRob3VnaC4gIFRoZSAqbWdjdiogcGFja2FnZSB3aWxsIHByb3ZpZGUgYSBiYXNpYyBwbG90IGFzIGZvbGxvd3M6CgpgYGB7ciBtZ2N2X3Bsb3R9CnBsb3QobW9kX2dhbTEpCmBgYAoKCldoYXQgeW91J3JlIHNlZWluZyBpcyBhICpjb21wb25lbnQgcGxvdCogd2l0aCB0aGUgeS9maXR0ZWQgdmFsdWVzIGNlbnRlcmVkLiAgQWxtb3N0IG5vIG9uZSBmaW5kcyB0aGlzIGludHVpdGl2ZSBpbiBwcmFjdGljZSwgbm9yIGlzIGl0IHVzdWFsbHkgd2hhdCB0aGV5J2QgcmVwb3J0LiAgSW5zdGVhZCB3ZSBjYW4gZ2V0IHByZWRpY3Rpb25zIGF0IGtleSB2YWx1ZXMgb2Ygb3RoZXIgY292YXJpYXRlcyAoZS5nLiBoZWxkIGF0IHRoZWlyIG1lYW4gb3IgcmVmZXJlbmNlIGNhdGVnb3J5KS4gIEluIHRoaXMgY2FzZSB3aXRoIG9ubHkgb25lIHByZWRpY3Rvciwgd2UganVzdCBnZXQgdGhlIHByZWRpY3RlZCB2YWx1ZXMgb24gdGhlIG9yaWdpbmFsIHNjYWxlIG9mIHRoZSBvdXRjb21lLiAgV2UgY2FuIHVzZSB0aGUgYGdnZWZmZWN0c2AgcGFja2FnZSBmb3IgdGhpcyAob3Igc2VlIG15IFt2aXNpYmx5IHBhY2thZ2VdKGh0dHBzOi8vbS1jbGFyay5naXRodWIuaW8vdmlzaWJseSkpLgoKYGBge3IgdmlzdWFsaXplX2luY29tZV9tYXJnaW5hbF9lZmZlY3R9CmxpYnJhcnkoZ2dlZmZlY3RzKQoKcGxvdF9kYXQgPC0gZ2dwcmVkaWN0KG1vZF9nYW0xLCB0ZXJtcyA9ICJJbmNvbWUiKQoKZ2dwbG90KHBsb3RfZGF0LCBhZXMoeCA9IHgsIHkgPSBwcmVkaWN0ZWQpKSArIAogIGdlb21fcmliYm9uKGFlcyh5bWluID0gY29uZi5sb3csIHltYXggPSBjb25mLmhpZ2gpLCBhbHBoYSA9IC4yNSkgKwogIGdlb21fbGluZShjb2xvciA9ICdkb2RnZXJibHVlJykgKyAKICBsYWJzKHggPSAnSW5jb21lJykKYGBgCgojIyMgTW9kZWwgQ29tcGFyaXNvbgoKV2UgY2FuIGNvbXBhcmUgbW9kZWxzIHZpYSBBSUMgb3IgdGhlIEdDViBjb2VmZmljaWVudCwgYm90aCBvZiB3aGljaCBmb2N1cyBvbiB0aGUgcHJlZGljdGl2ZSBjYXBhYmlsaXRpZXMgb2YgdGhlIG1vZGVsIG9uIG5ldyBkYXRhLgoKYGBge3IgbW9kZWxfY29tcGFyaXNvbn0KQUlDKG1vZF9sbSwgbW9kX2dhbTEpCmBgYAoKTGlrZWxpaG9vZCByYXRpbyB0ZXN0IChhcHByb3hpbWF0ZSkuCgpgYGB7ciBhbm92YV9nYW19CmFub3ZhKG1vZF9sbSwgbW9kX2dhbTEsIHRlc3Q9IkNoaXNxIikKYGBgCgoKIyMgTXVsdGlwbGUgUHJlZGljdG9ycwoKCiMjIyBMaW5lYXIgRkl0CgpgYGB7ciBtb2RfbG0yfQptb2RfbG0yIDwtIGdhbShPdmVyYWxsIH4gSW5jb21lICsgRWR1ICsgSGVhbHRoLCBkYXRhPXBpc2EpCnN1bW1hcnkobW9kX2xtMikKYGBgCgojIyMgR0FNCgpgYGB7ciBtb2RfZ2FtMn0KbW9kX2dhbTIgPC0gZ2FtKE92ZXJhbGwgfiBzKEluY29tZSkgKyBzKEVkdSkgKyBzKEhlYWx0aCksIGRhdGE9cGlzYSkKc3VtbWFyeShtb2RfZ2FtMikKYGBgCgpXZSBjYW4gc2VlIHRoZSBmaXQgYXBwZWFycyB0byBiZSBiZXR0ZXIuICBJbiBhZGRpdGlvbiwgYXMgdGhlIGVmZmVjdGl2ZSBkZWdyZWVzIG9mIGZyZWVkb20gaXMgMSBmb3IgSGVhbHRoLCBpdCBpcyBlc3NlbnRpYWxseSBhIGxpbmVhciBmaXQsIHNvIGlmIGRlc2lyZWQsIHdlIGNhbiBsZWF2ZSBpdCBhbW9uZyB0aGUgYHBhcmFtZXRyaWNgIHRlcm1zLgoKYGBge3IgbW9kX2dhbTJfcGxvdH0KIyBwbG90KG1vZF9nYW0yKQoKbGlicmFyeShwYXRjaHdvcmspICMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJ0aG9tYXNwODUvcGF0Y2h3b3JrIikKCmcxID0gCiAgZ2dwcmVkaWN0KG1vZF9nYW0yLCB0ZXJtcyA9ICJJbmNvbWUiKSAlPiUgCiAgZ2dwbG90KGFlcyh4ID0geCwgeSA9IHByZWRpY3RlZCkpICsgCiAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSBjb25mLmxvdywgeW1heCA9IGNvbmYuaGlnaCksIGFscGhhID0gLjI1KSArCiAgZ2VvbV9saW5lKGNvbG9yID0gJ2RvZGdlcmJsdWUnKSArIAogIGxhYnMoeCA9ICdJbmNvbWUnKQpnMiA9IAogIGdncHJlZGljdChtb2RfZ2FtMiwgdGVybXMgPSAiRWR1IikgJT4lIAogIGdncGxvdChhZXMoeCA9IHgsIHkgPSBwcmVkaWN0ZWQpKSArIAogIGdlb21fcmliYm9uKGFlcyh5bWluID0gY29uZi5sb3csIHltYXggPSBjb25mLmhpZ2gpLCBhbHBoYSA9IC4yNSkgKwogIGdlb21fbGluZShjb2xvciA9ICdkb2RnZXJibHVlJykgKyAKICBsYWJzKHggPSAnRWR1JykKZzMgPSAKICBnZ3ByZWRpY3QobW9kX2dhbTIsIHRlcm1zID0gIkhlYWx0aCIpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0gcHJlZGljdGVkKSkgKyAKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IGNvbmYubG93LCB5bWF4ID0gY29uZi5oaWdoKSwgYWxwaGEgPSAuMjUpICsKICBnZW9tX2xpbmUoY29sb3IgPSAnZG9kZ2VyYmx1ZScpICsgCiAgbGFicyh4ID0gJ0hlYWx0aCcpCgooZzIgKyBnMyArIGcxICsgcGxvdF9sYXlvdXQobnJvdyA9IDIpKSAqIHRoZW1lX3RydWVNaW5pbWFsKCkKYGBgCgoKIyMjIDJEIFNtb290aHMKCmBgYHtyIG1vZF9nYW0zfQptb2RfZ2FtMyA8LSBnYW0oT3ZlcmFsbCB+IEhlYWx0aCArIHRlKEluY29tZSwgRWR1KSwgZGF0YT1waXNhKQpzdW1tYXJ5KG1vZF9nYW0zKQpgYGAKCgpgYGB7ciBtb2RfZ2FtM19wbG90fQojIHJlcXVpcmVzIHRoZSB2aXNpYmx5IHBhY2thZ2UsIG90aGVyd2lzZSBzZWUgdGhlIG1haW4gZG9jdW1lbnQ7IAojIGRldnRvb2xzOjppbnN0YWxsX2dpdGh1YigibS1jbGFyay92aXNpYmx5IikKdmlzaWJseTo6cGxvdF9nYW1fM2QobW9kZWwgPSBtb2RfZ2FtMywgbWFpbl92YXIgPSBJbmNvbWUsIHNlY29uZF92YXIgPSBFZHUsIHBhbGV0dGU9J2JpbGJhbycsIGRpcmVjdGlvbiA9IDEpCmBgYAoKIyMjIE1vZGVsIENvbXBhcmlzb24KCmBgYHtyIG1vZGVsX2NvbXBhcmlzb25fcmVkdXh9CkFJQyhtb2RfbG0yLCBtb2RfZ2FtMikKYGBgCgo=